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Abstract 

Homogenization of the equations of motion for a three dimensional periodic elastic 
system is considered. Expressions are obtained for the fully dynamic effective material 
parameters governing the spatially averaged fields by using the plane wave expansion 
(PWE) method. The effective equations are of Willis form jT] with coupling between 
momentum and stress and tensorial inertia. The formulation demonstrates that the 
Willis equations of clastodynamics are closed under homogenization. The effective 
material parameters are obtained for arbitrary frequency and wavenumber combina- 
tions, including but not restricted to Bloch wave branches for wave propagation in 
the periodic medium. Numerical examples for a ID system illustrate the frequency 
dependence of the parameters on Bloch wave branches and provide a comparison with 
an alternative dynamic effective medium theory [2] which also reduces to Willis form 
but with different effective moduli. 

1 Introduction 

The effective medium concept has proven to be very useful for modelling composite ma- 
terials with periodic microstructure. Efficient methods exist for computing the effective 
elastic moduli under static and quasistatic conditions. A more challenging task is to define 
frequency-dependent dynamic effective constants which are capable of describing phononic 
band gaps of the wave spectrum [3] at the expense of non-classical features such as negative 
[1] and tensorial [5] density. Effective material models that can retain a dispersive signature 
of the initial material microarchitecture have significant potential application in the design 
of sonic metamaterials. The purpose of this paper is to provide an analytical formulation 
for the effective material parameters of a periodic elastic medium at finite frequency. 
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A general framework for describing dynamic effective medium theories has been devel- 
oped by Willis, beginning with variational arguments for non-local behavior in mechanical 
[H [7] and electromagnetic systems [8]. The predicted constitutive relations modify both 
the elasticity [5] and the inertia [TU] , and include coupling between stress / strain on the one 
hand and momentum/velocity on the other, see pQ for a review. The notion of anisotropic 
inertia has appeared in several mechanical contexts, such as stratified layers of fluids [llj 
and elastic composites with strongly contrasting constituents [12] . Coupling between parti- 
cle velocity and stress, although not usually considered in constitutive theories for continua, 
is to be expected for theories relating spatial means of these quantities because of the under- 
lying inhomogeneity of the displacement and stress fields within the region being averaged. 
Interest in the Willis constitutive model has intensified with the observation [5] that the 
non-classical material properties could be realized, in theory at least, by suitably designed 
microstructures with internal springs, masses, and gyroscopes. Derivations of the Willis 
equations have been demonstrated for periodic systems in ID [T21 [21 U3] and 3D [T2], and 
the structure of the equations has been rigorously proved for both electromagnetic and 
elastic materials [16j . 

Despite the growing awareness that the Willis constitutive model is the proper dynamic 
effective medium formulation for periodic systems there is as yet no simple means of cal- 
culating the effective material parameters at finite frequency. Several different techniques 
have been proposed to address this deficiency, none of which is easily implemented for 
arbitrary 3D problems. Thus, closed form expressions were obtained in [13] for ID strati- 
fication using Floquet modes, where the latter have to be determined numerically. A more 
general numerical technique has recently been developed and applied to ID [TJ] and 3D 
|15] composites. This approach employs a background or reference material, and solves 
for the polarization field (difference in stress/momentum between the actual and that of 
a comparison body [17]) to arrive at expressions for the effective parameters in terms of 
system matrices. The numerical procedure, which uses a combination of plane wave expan- 
sions and discretization schemes, is complicated by the required selection of a comparison 
material, the choice of which has no bearing on the final outcome. A quite different tech- 
nique [2] suited to lD-periodically stratified elastic solids uses the monodromy matrix (i.e., 
propagator for a single period) to explicitly define effective material parameters from its 
logarithm. All of these methods derive the density, stiffness and the coupling parameter, S, 
as functions of frequency and wavenumber. This means, perhaps surprisingly, that equa- 
tions of Willis form but with different parameters can be obtained for the same medium, 
depending on the homogenization scheme employed, an aspect discussed in £}6j 

Our purpose here is to develop, in so far as possible, analytical expressions for the 
dynamic effective material parameters of a periodic elastic material. The homogenization 
scheme seeks the constitutive relations and the equations of dynamic equilibrium governing 
the spatial average of the periodic part of the Bloch wave solution. The initial heterogeneous 
system is here assumed to be defined by pointwise constitutive equations of the form 
proposed by Willis, which includes 'classical' elastodynamics as the case of isotropic density 
and zero coupling between stress and momentum. We show that the homogenized equations 
are also of Willis form with dynamic effective material parameters that are non-local in 
space and time, thus demonstrating that the Willis constitutive theory is self-consistent 
and closed. The homogenized effective parameters are conveniently expressed in terms 
of plane wave expansions (PWE) of the original material parameters, and accordingly we 
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call this the PWE method. This approach allows us to define effective constants for any 
frequency-wavenumber combination, including, but not restricted to, values of {u, k} on 
the Bloch wave branches. 

We begin in <J5] with an overview of the Willis constitutive equations followed by a 
summary of the PWE homogenization results. Section Sj3] is devoted to the derivation of 
the PWE effective material parameters for a scalar wave model that is simpler than the 
fully elastic problem, but which exhibits the main features of the homogenization scheme. 
The elastodynamic effective parameters are obtained in $3] and some of their properties 
are examined in £}5j Examples of the PWE effective moduli are provided in £|6] based on 
numerical calculations for a ID system, and comparisons are made with the parameters 
predicted by the monodromy matrix (MM) homogenization scheme [2]. Conclusions and 
final words are given in £J71 

2 Problem setup and summary of the solution 

2.1 Willis elastodynamic equations for a heterogeneous medium 

The field variables are the vectors of displacement, momentum and body force, u, p, and 
f; the symmetric second order tensors of stress, strain, and strain source, <r, £, 7, with 
£ = ^[Vu + (Vu) T ]. They are related by the system of elastodynamic equations which 
we take in the Willis form (in order to demonstrate in the end that this model is closed 
under homogenization for periodic materials). Accordingly, the material parameters are 
the second order tensor of mass density p, the fourth order elastic stiffness tensor C, and 
a third order coupling tensor S. The components of the above material tensors in an 
orthogonal basis of the coordinate space satisfy the symmetries 

Pij = Pjii Cijki = C^nj, Cijki = Cjiki, Sijk = Sjik- (1) 

The two former identities imply that the density and stiffness tensors are Hermitian: p = p + 
and C = C + (* and + hereafter denote complex and Hermitian conjugation). Assuming 
both p and C positive definite ensures positive energy ([3]), see below. Let S be generally 
complex. Field variables are functions of position x and time t, while the material tensors in 
the initial heterogeneous medium are considered as functions of x only (the homogenization 
leads to the PWE effective parameters which, after inverse Fourier transform, are functions 
of both x and t). 

The equation of equilibrium and the generalized constitutive equations in Willis form 
[3 [TU] are, respectively, 



p = diver + f , 




where the dot implies time derivative. Note that the terms of ([2]) with S taken in com- 
ponent form imply <r»j 3 S^j) k u k and pi 3 S^ jk ^ (e - j) jk = -Sy^ (e - j) jk , where the 
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permutable indices are enclosed in parentheses (see §5. II for further details). By ([I]) and 
([2]) with 7 = 0, the time-averaged energy density is 

E = ^Re(<r + e + u+p) = - (e + Ce + u + pu) . (3) 

It is obvious that eqs. ([2]) include 'classical' linearly anisotropic elasticity with pij = p5ij, 
Cijki = Ckiij = Cjiki and Siji = 0. 

Equations (J2J) may be solved to find u, p, and a for given forcing functions f and j. 
Accounting for the effect of two driving-force terms is an important ingredient of homoge- 
nization concepts [16]. In particular, dependence of the solutions on the forcing terms will 
turn out to be crucial in finding the PWE homogenized equations. 

2.2 The periodic medium 

We consider the material with T-periodic parameters: 

M x +Ej=i™; a j) = M x ) for h = c, P , s, (4) 

where x € M. d (d =1, 2 or 3), rij £ Z, and the linear independent translation vectors aj € M. d 
define the irreducible unit cell T = X^/=i *i a i e IPj 1]) °f ^ ne periodic lattice. Let {ej} 
with j = 1, . . . , d be an orthonormal base in M d . Denote 

By = Aej, by = (A~ 1 ) T e j , r = {g : g = Y'J {1-iijbj. Uj G Z}, (5) 

where • = 5jfc (• is the scalar product in W 1 ), T means transpose, and g = X^=i ffi e j 
is an element of the set T of reciprocal lattice vectors whose components {gj} in the base 
{ej} take all values 27r(A~ 1 ) T Z d . The Fourier transform maps a T-periodic function 
h(x) to a vector h in the infinite-dimensional space Vg associated with g £ T, where the 
components /i(g) of h are Fourier coefficients of h(x) : 

h(g) = IT]" 1 r /i(x)e-*' x dx ^ /»(x) = E ge r^(g)e ig x (6) 

JO 

(all quantities in the Fourier domain are indicated by a hat hereafter). In particular, the 
average over the single cell is defined by 

(h) = IT]" 1 f /i(x)dx ( = h(0) = e+h where e (g) = <5 g0 ). (7) 
Jo 

Practical calculations in the Fourier domain imply truncation of vectors and matrices in 
Vg to a finite size, i.e. reducing V s to finite dimension. This is tacitly understood in 
the following when using the concepts of matrix determinant and inverse as convenient 
shortcuts. 
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The governing equations (|2|) with periodic coefficients ([!]) will be solved for both the 
wave field variables and the force terms in Bloch form 

/i(x,t) = /i(x)e i(kx ^' ) , h(x) T-periodic, for h = u,p,o-,£,f,7, (8) 

where /i(x) implies the periodic part of the full dependence of h on x (this convention 
reduces subsequent notation), and the frequency u and wave- vector k are assumed real- 
valued. Note that one may think of f and 7 in the form ([8]) as phased forcing terms which 
'drive' the solution. 



2.3 Summary of PWE homogenization results 

For h expressed in Bloch wave form, (|8|), define the effective field variable 

/i cff (x,t) = (h) e^-^ . (9) 

Using the effective sources (f eff ,7 eff ) = ((f), (7}) e H k ' x_a '*) implies ignoring the influence 
from the forcing f(x), 7(x) with zero average over a single cell, which is natural for the 
homogenization theory aimed at recovering the effective wave fields of the form @ . One of 
the principal results of the paper is that, given the periodic medium with the Willis form 
([2]) of the governing equations (which incorporates 'classical' elastodynamics model), the 
equations describing relations between the effective wave fields u eff , c eff , p eff and e eff = 
^[Vu cff + (Vtt cff ) T ] and the effective forcing f cff and j cS are also of Willis form: 

C cff gcff\ / c ff _ 7 cffN 

p cff = diva eS + f cff , [ 

„ S cff+ pCS J 




The effective parameters C , p and S are non-local in both space and time (despite 
the fact that the exact parameters of the original periodic material are stationary). Corre- 
spondingly, they are functions of frequency u) and wavenumber vector k in the transform 
domain (g§-, — > ikj, —ioo), where their explicit expressions are as follows: 

qf w (w, k) = (C ijkl ) - q^G^qfc^, (Ha) 

Sg k (u>, k) = (S ijk ) - q^G^V, (lib) 

Pik ( w > k ) = (Pik) + r&Gpqiqk, (He) 

with the summation on repeated indices hereafter implicitly understood. (Noting that the 
indices of tensors C and S are appropriately split between C 3 B u and M. d 3 x, we assume 
for brevity that d = 3 and so all roman indices run from 1 to 3). The tensors q, f and 
Gr\° are defined on V^o®C 3 where is the Fourier vector space restricted to g € T\0. 
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Each 'element' q^-fc, and &° q with fixed indices i...q is, respectively, a vector and a 
matrix in V^o- Their components are as follows: 

^(g) = (ffl + ^i(g)-4fcfe)' ger\o ; (lid) 

fik(s) = upik(g) + (9j + kj)Sijk(g), g € r\0; (lie) 
J] g/ Gj [g, g'] [g', g"] = 5 lk 5 gs „ with g, g', g"e T\0, and (llf) 

z-l [g, g'] = (s-i + + k)C ijk i(g - g') - ^ 2 A.fc(g - g') 

- + kj)S ijk (g - g') - utfj + %)5^(g' - g), (llg) 

where Sijk = S^j^ = Sjik- Note that since the matrix is Hermitian (G^[g, g'] = 

Gp°* [g',g]), the effective material tensors (fTTj) established in this paper recover the sym- 
metries ([T]) of the corresponding tensors of the initial periodic material (although not the 
sign definiteness, see £15.11 for more details). The Willis formulation of elastodynamics 
equations is therefore closed under the PWE homogenization. 



3 PWE effective parameters: derivation for a scalar wave 
model 

Instead of proceeding with the elastodynamic equations ([2]) in full, we first provide a 
detailed derivation for the case of a scalar wave equation. The scalar case is more revealing 
as it avoids the multitudinous suffices of the equations of elasticity while retaining the 
essence of the homogenization. We therefore ignore the obvious fact that an uncoupled 
acoustic scalar wave in solids with x-dependent material parameters is restricted to x € M. d 
with d =1 or 2 since the scalar problem is purely a methodological shortcut. Once the 
scalar derivation has been completed, it is straightforward to generalize the results to the 
original elastodynamics system ([2]) in x G 1R , d =1, 2 or 3, see $H 

3.1 Exact solution of the scalar wave problem 

For the scalar wave model, the displacement u, momentum p, body force / and density 
p are scalars, while the stress a, strain e = Vn, strain source 7, and coupling parameter 
S are vectors and the 'shear' stiffness /x is a d x d matrix in C d (disregarding the trivial 
case of d = 1). According to (pQ), p is real and fj, Hermitian while S may be complex. Let 
these material parameters be T-periodic in the sense (j4|) and take the wave variables and 
sources in Bloch form (JH]). Then the scalar- wave version of the governing equations (|2|) in 
component form, with j = d/dxj, is 

/o-A / fiji SA /uj + ikiu - 7A 
-iup = f + ik j a j + a jJ , 1= J J. (12) 
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Applying Fourier transform Q reduces (fT2|) to the following algebraic equations in V g : 

-iwp(g) = /(g) + i(gj + k^ajig), 

( Hlis ~ g') ^(g - g')\ Kh + g'Mg') - 7,(g')\ (13) 

§' er Wr(g-g') p(g-g')/ \ -^n(g') / 

It is convenient to pass to a matrix form in Vg while keeping explicit the coordinate indices 
in x-space (this makes the eventual generalization to the fully elastodynamic problem 
clearer). Recall the notation h for the vector in Vg with the components /i(g) (see ([6])), 
and introduce the matrices M_,7, D/,N, V; (g Vg <g> V s when j,l are fixed) with the 
components 

Mji [g,g'] = Aji(g - g'). A [&g'] = + **)<W, (14) 
# [g,g'] = p(g - g'), Vi [g,g'] = 5 z (g - g')- 

Note that Mji[g',g] = My [g, g'] etc., i.e. the matrices Mjj, D; and N are Hermitian, 

while V; is not once Si (x) is complex. Equations (|13p can now be expressed in a matrix 
form on V s as 



-iup = f + iDj&j, 




-V+ N / \ -iuu 




(15) 



where = M j7 D ; - u)V j} R = V+D^ + wN. 



Reduction to a single equation for the displacement provides a definition of the impedance 
matrix Z in the Fourier domain: 

Z u = f — where 

\ 3 , , , „ (16) 

Z(w, k) = DjMjiDj - w(VtDj- + Dj-Vj) - w 2 N ( = D^Qj - uR) = Z + (w, k). 

We note that the condition for existence of free waves in the absence of sources, det Z(w, k) = 
0, defines the Bloch eigen- frequencies u) (k). Our interest is ultimately in homogenized equa- 
tions of motion allowing for the Bloch waves. It is however simpler to first proceed under 
the assumption that Z(w,k) is not singular (det Z ^ 0). The issue of how to handle a 
singular Z is resolved afterwards in §3.31 by introduction of a slightly modified impedance. 

Assuming Z is non-singular the Green's matrix function G(w, k) is defined as the inverse 
of Z. The exact scalar solution then follows from (1161) as 



u = Gf - iGQ + 7i , with G = Z" 1 ( = G+) . (17) 
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3.2 Homogenized equations 

Taking the average ([7|) of the wave solution (fT7|) and the constitutive equations (fl~5|) 9 yields 

(it) =e + Gf-ie + GQ+7 j , 

(a,) = ie+Q.-Gf + e+(Q i GQ+ - M jt )% (18) 
(p) = -*e + RGf - e+(RGQ+ - Vf)% 



Our objective is to rearrange eqs. (|18p into a form relating the effective wave fields to the 
effective forcing terms with the latter viewed as independent variables. As a first step in 
this direction, we insert the effective sources (see eq. ([9])) 7° ff = (7j)e and f cff = (/) e in 
eqs. (|18|) and rewrite the result in a compact form 

(u)=a(f) -iB* ( 7j ), 

(a j ) = iB j (f)-A jt < 7z ), (19) 
(p) = -ib* (f) +a* 

where the scalars a, b € C, the vectors a, @ € C d and the matrix A € C d (8> C d are 

a = e + Ge, b = e+GR+e, 8- = e + Q 7 Ge, 

. . . 3 3 \ , . (20) 

Aji = (fiji) - e+QjGQ+e, aj = (Sj) - e+Q.GR+e. 

Note that a is real and A + = A. Eliminating (/) using (]19P i allows us to recast eq. ([191b 
and (|19p q in a form reminiscent of the Willis constitutive relations, 






(21) 



where Xj and Y are at this stage arbitrary. Comparing (12ip with the assumed structure 
of the effective constitutive relations applied to fields in the presence of a strain source, eq. 
(IIOP ? and bearing in mind ([9|) yields 



(22) 





4 

-sf 

Equating the contributions to (aj) and {p) in (|2ip and ()22[) from 7 gives, in turn, 

M |f = A iZ + i/3,/3f , Sf = a, + h -B r (23) 

Full equivalence requires Y = p cff and X = S cff , which implies two more identities found 
by equating the contributions from (u) to (p) and to (cry), respectively: 

ujp cS =a~ 1 b* -Sf*kj, ojSf = fifki-a^Bj. (24) 
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The first equation provides an expression for p and the second is an alternative to ([23]) 9 
for S eff (equivalence of ([25]) 9 and (f24"]) 9 is noted below). 

Combining eqs. (|23p and (|24p i and reinstating the notations from (|20p yields the effec- 
tive parameters in the form 

ff/ , v , . . . , (e + Q,Ge)(e+GQ+e)+ , 

Hf(u, k) = - e+Q,GQ+e + 1 - ^ - - , (25a) 

e+Ge 

* « « Ip+Gr+pI 2 

p eff (w, k) = <p) + e+RGR+e - l - — L (25b) 

e+Ge 



(e+GR+e)(e+GQ+e) 

Sf(u, k) = (Sj) - e+Q.GR+e + ^4 i. (25c) 

e+Ge 



In the derivation of ([25]) . we have taken note that Z = DjQj — wR = Z+ and used identities 
such as 

k • ft* = fye+GQje = e+GQ+D,e = e+G (z + wR + ) e = 1 + we+GR+e. (26) 

The equivalence of the expressions (fjSlb and (fjMh for S eff follows from using the explicit 
formulas (]25j) and identities of the above type. 



3.3 Regularized form of the effective parameters 

By definition (|17p . the Green's function G(cj,k) = Z _1 (o;,k) diverges on the Bloch 
branches Qb (3-parameter manifold in {oj, k}-space) defined by the Bloch dispersion equa- 
tion 

B 9(w,k) B : detZ = O Zu B = 0, (27) 

where the null vector ub of Z is the Bloch eigen-mode in the Fourier domain. At the same 
time, the effective dynamic parameters given by eqs. ()25[) generally remain finite on f^B 
(accidental exceptions such as degenerate and certain other points on Qb are disregarded 
for the moment and will be discussed in £ j3.4p . This can be shown by inspecting the limit 
of eqs. (I25p as lj, k tend to a (single) point (lj, k) B on a Bloch branch, from which it is 
seen that the right-hand side members of (I25p that diverge on f^B in fact compensate and 
cancel each other. However, such limiting behavior of eqs. (]25p prevents their numerical 
implementation by evaluating each member independently before summing them up. The 
difficulty can be circumvented if we recast eqs. (|25p in the analytically equivalent but 
explicitly different form that is rid of the terms diverging on Qb- Since the divergence in 
(|25p is due to Green's function G, the idea is to extract the part G reg of G which is regular 

on f^B- Consider such a partitioning of G in the general form G (w, k) = G reg + e^vv" 1 ", 
where v — > ub and e — > as u), k — >• (w,k) B . It is not unique. For instance, one standard 
way is to invoke the spectral decomposition 

Z (lj, k) = Z' + ?vv+, G (w, k) = G' + ^ 1 vv+, (28) 
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where ? and v are the corresponding eigenvalue and eigenvector of Z (= Z + ) which satisfy 
v — > ub and ? — > as c*>, k — )• (w, k) B , and G' = Z' (the pseudoinverse of Z) is regular 
on Ob- Substituting (pHjb for G in (f25j) removes the diverging terms and yields an explicit 
form expressed through G'. Unfortunately, explicit calculation of G' is relatively laborious. 
In this regard, we advocate the use of another partitioning of G into regular and diverging 
parts which requires merely singling out Fourier components with zero and non-zero g. 
This is an essential ingredient of our homogenization method that ensures both robust and 
straightforward numerical implementation. 

With the above idea in mind, denote by V s ^o the vector space that is associated with 

the set T\0, i.e. is orthogonal to the vector e (see ([!]) and (j7|)). Let Z\ be the restriction 
of Z over this space V s ^o and G\o be the inverse of Z\ on V s ^o : 



z \o 



^ Q ; G\ = Z\~o ( for det Z\o + 0) (29) 



(note that G\o as defined is not a restriction of G = Z 1 on V^o)- Introduce the vectors 
w, q j; f G Vg^o, 

w = kjC[j — ur where 

q, = Q+e - (eQ+e)e, q 3 (g) = ( 9 i + fc,)£* (g) " ^* (g) , for (30) 
r = R+e-(e+R+e)e, r (g) = w£(g) + ( 5j + (g) , 

Also denote (Z) = (Ze) and (G) = (Ge) which is a slight generalization of the average © 
to include the projection of a matrix onto g, g/ = 0. Thus 



(31) 



(Z) ee Z [0, 0] = e+Ze = kjki(nji) - lo z { P ) - ukj{Sj + S*), 
(G) ee G [0, 0] = e+Ge = det Z\„ det Z _1 . 
Using the notations (j30|) and ()3ip . we note the identical decompositions 

Z(w,k) = Z\ + we + + ew + + (Z)ee + , det Z = (Z) det Z\ - w + Z\ w, (32) 
where Z\ is adjugate of Z\ . Hence if det Z\ ^ as assumed in (j29]) 9. then 

{G)- 1 = (Z) - w+G\ w, (33a) 

det Z (u, k) = (G)- 1 det Z\„ (= ^ (G)" 1 = on fi B ) , (33b) 

Zv = (G)" 1 ^ with v (w, k) = e — G\ w (— > ub asw,k-> S7b) j (33c) 

G(w,k) = Gy + (G)vv+, (33d) 

where eq. (]33dp yields the sought-for partitioning G into two parts, one regular and the 
other diverging on Bloch branches Qb- Note that generally (i.e. for G\qW 7^ 0) the quan- 
tities G\ and v (— > ub) in (|33d]) are different from, respectively, the pseudoinverse G' 



and the eigenvector v (— > ub) in (|28D . 
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Writing (|33dp in the transparent form 

G\ = G — (e + Ge) _1 Gee + G (34) 

and using the notations (f30|) . reduces eqs. (f25j) to the simpler form 

Mif(w,k) = (/ij-j) - q+G\ q,, (35a) 
Sf(w,k) = {S j )-qfG\ f, (35b) 
p eff ( W ,k) = (p)+f + G\ f. (35c) 

It is emphasized that the effective parameters explicitly defined in (|35p are by construction 
the same as those in (|25|) . At the same time the terms in (|25|) divergent on f?B are 
eliminated in (|35p . in this sense (|35[) are uniformly convergent. The new expressions (|35p 
are determined from the regular part G\o of the Green's function associated with g G T\0 
which is straightforward to compute; moreover, their explicit form is simpler than that of 



3.4 Exceptional cases 

The expressions for the effective parameters (|35p defined by the matrix G\ = Z^ 1 are 

obtained for the general case of non-singular Z\ - If Z\ happens to be singular, then 
numerical implementation of (|35p produces diverging terms (though eqs. (|25p expressed 
via G = Z remain finite for detZ\ = so long as detZ ^ 0). The properties of the 
wave solutions that occur specifically at values of cj,k rendering Z\ (w,k) singular are of 
interest. Denote the set of these exceptional points (3-parameter manifold in {u, k}-space) 
by Q 

exc so that 

Q cxc 3 (w,k) cxc : detZ\ o = 0. (36) 

Note that eqs. (|33p are invalid on O exc . 

Consider first 'generic' points of 0, cxc which lie off the Bloch branches (|27|) . i.e. (w, k) cxc ^ 

B (=> detZ / 0). Then by pi^ o (G) = (Ge) = and hence the averaged wave (fl^D i 

is (u) = 'i(Z) -1 7jW + Z\oqj- It is generally nonzero, which is not surprising as (u) = 
would correspond to the zero effective response to the non-zero effective forcing. At the 
same time, the wave response to a pure force source (7^ = 0) at the points (w,k) cxc ^ 
has a zero mean over the unit cell T in x-space. It is therefore somewhat satisfying and 
physically sensible that the effective parameters (j35|) indicate such exceptional points by 
becoming numerically divergent. 

Now assume a generally non-empty set £l cxc H J^b (2-parameter manifold in {i(j,k}- 
space) of the points of Bloch branches (j27|) at which Z\ happens to be singular. By (f32j) 9. 

simultaneous equalities det Z = and det Z\ = yield w + Z\ w = 0. It is natural to 

assume that Z\ w is nonzero. Then it is a null vector of Z\ and it is orthogonal to w. 

Hence, by (f32l) i . Z\ w is also a null vector of Z, i.e. it is a Bloch mode ub- Thus, if 
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det Z\ = incidentally occurs at some point of a Bloch branch, then the Bloch mode at 

this point is ub = Z\ w, which is by construction orthogonal to e and so its original ub (x) 
in x-space has zero mean (ub) = over the unit cell T. Let us call it a zero-mean Bloch 
mode. This is in contrast to 'normal' points of SIb where the Bloch modes are defined by 
(|33c|) as ub = e — G\ w and thus have non-zero mean (ub) = 1- Since zero- mean Bloch 
modes are zero (trivial) solutions of the effective equations for free waves, their explicit 
exclusion from the framework of the homogenization theory, as signalled by divergence of 
the effective material parameters (|35p . is physically consistent. Note that the divergence 
at the points of zero- mean Bloch modes is 'genuine' (analytical), i.e. it occurs regardless 
of which explicit expressions are used for the effective parameters (unless more restrictions 

are imposed apart from r2 exc H ^b, such as Z\ w = 0). 

In conclusion, we note another exceptional case for the effective parameters (|35|) is the 
possible existence of double points (self-intersections) on Bloch branches Qb- It can be 
shown, using e.g. the spectral decomposition and the pseudoinverse of Z, that a double 
point on Qb causes a genuine (in the sense mentioned above) divergence of the effective 
parameters, which conforms to the fact that a double point admits the construction of a 
zero- mean Bloch mode (ub) = from a pair of null vectors of Z. 



3.5 Summary for a scalar wave problem 

The equations governing the averaged effective fields u cS , a cS , p eS and £ cff = Vu cff and 
their relation to the averaged applied body force f cS and strain source 7 (see ([!])) are of 
Willis form: 



p cS =diw cff +/ cff , 



eff gcffx / off .yrfF 



S cff+ p cff 



n cff 



(37) 



where p eS , uff and Sj are non-local in space and time. They are defined in the transform 
domain by eqs. (|35p with the following explicit form 



Sf (w,k) 



> = < 



{Hi) 
(Sj) 
(P) 



+ g , g ?r\o 6 \o [S'S'] X< 



fijn(-g)(k n + 9n) ~ uSjis) h P (s')(k p + g' p ) - uS t (g! 
Hl(-g)(k + 9i) - w5j(-g)] [up(s') + (K + 9n)S n (s') 
wp(-g) + (kj + g^Sji-g)] Lp(g') + (h + gfiSiig! 



(38a) 
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where G\ = Gt is the inverse to the matrix Z\ with the components 

Z\o [g, g'] = Ukj + + gDfijiig - g') - w 2 /5(g - g') 

-wfo + fc)4-(g - g') - wfo + ^)S*(g' - g)) . (38b) 

Note that p eS = /i eff+ , /9 cff is real and Sj S is complex- valued, see ^5.11 

The effective parameters defined by eqs. ([38]) are regular at any uj, k excluding possible 
double points (intersections) of Bloch branches, and the exceptional points (w, k) cxc G O ex c 

where det Z\ = (a numerical divergence appears in the vicinity of such points). If a point 
(w,k) exc incidentally occurs on a Bloch branch then the corresponding Bloch eigen-mode 
ub (x) has a zero mean (ub) = over the unit cell T. Zero- mean Bloch mode can also 
exist at the double point of a Bloch branch. Barring these exceptions, eqs. (|38p are regular 
on Bloch branches which are defined by the dispersion equation (see (j33a|) ) 

(G)- 1 (w,k) = (detZ\ o ^0). (39) 

The term (G 1 )" 1 = Z cS which defines the dispersion relation is expressed as a more physi- 
cally meaningful quantity in §5.21 



4 Elastodynamic effective parameters 

The procedure of adapting the scalar acoustic model to the vector elastodynamic model 
requires that the scalar and vector variables u, p, f and c, £, 7 of the acoustic case become, 
respectively, vectors and symmetric second order tensors in C 3 . Accordingly, the scalar 
density p, coupling vector S and the stiffness matrix fj, increase by two orders as tensor 
and become p, S and C defined in H2.ll The component form of the various quantities for 
the two problems are related as follows: 

scalar system elastodynamic problem 

U, P, f, (Tj, £j, Jj, p, Sj, pjl -4 Ui, Pi, fi, (Tij, Eij, Jij, Pik, Sijk, Cijkl- (40) 

Based on this equivalence the generalization of the scalar eqs. f)35 [) proceeds by first re- 
placing the single suffix j in the former with ij in the latter, then using definitions such as 
those in eq. (|30p to assign additional suffices, ultimately yielding the elastodynamic result 

Expressing the effective elastodynamic parameters (|11|) via the matrix ensures that 
they are generally regular on the Bloch branches. This is an important feature that deserves 
further scrutiny. By analogy with (|27p . denote the manifold of Bloch branches for vector 
waves by Qb (bold- lettered to distinguish it from Qb of the scalar case), so that 



fi B 9 (w, k) B : det Z = 44> Zu B = 0, 



(41) 
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where the matrix Z and the vector ub have components [g, g'] and uj^ (g') in V s x C 3 . 
If detZ X ° ^ and so G X ° = Z\ 0_1 is well-defined (see ([II])), then 

det Z (u, k) = det(G)" 1 det Z X °, (G)" 1 = (Z) - w+G X °w with det(G)" 1 = on n B - 

(42) 

This is similar to the scalar version (j33a[) - (|33bj) with the important distinction that now 
(Z) = Z[0,0] = e+Ze and (G) = G[0,0] = e+Ge are 3 x 3 matrices in C 3 and w 6 
C 3 <S> ^g^Oj e € C 3 <8> V s= q. Their components in C 3 are defined by 

(Z) : (Z) ik = kjki{Cij k i) - u 2 {p ik ) ~ ujkj{Si jk + S kji ), 

(G)" 1 = Z cff : Z* = (Z} lk - E gig , er \ < (8) (ft «0 *<r* («0 > ( 43 ) 
w, e : w ik = kjq ijk - u)r ik , e ik = 5 ik S g0 

(the quantity Z cff is expressed in more physical terms in §5.2|) . Note also that the vector- 
case generalization of ()33cp implies that the null vector of (G)^ 1 (= Z eff ) on fie is equal to 
the average (ub) of the Bloch eigenmode defined by (j4ip ?. Finally, consider the manifold 

fi e xc of exceptional points where det Z X ° = (cf. (j36|) ) and the explicit expressions (fTT|) 
diverge. For the points in Cl cxc f] CIb, i-e. lying off the Bloch branches, we have det(G) = 
(cf. (|42p i ) and hence the effective force f eff (x,i) Q with (f) € C 3 being a null vector of 
(G) produces a zero effective response (u) = 0. This is again similar to the scalar case 
(except that the vector (f) has to be specialized while (u) = V (/))■ If the intersection 
fiexc H Ob occurs, i.e. det Z^ = occurs on a Bloch branch, then, contrary to the scalar 
case, the vector Bloch mode does not generally have a zero mean (ub) = unless the null 

vector h of Z X ° also satisfies w + h = (see eq. (f3"2"j) i ). 

5 Discussion of the PWE effective model 
5.1 Properties of the effective constants 

As it was noted in £12.31 the effective material tensors given by eqs. (|lip retain the sym- 
metries (HJ) of the corresponding tensors of the initial periodic material. At the same time, 
the effective density and elasticity tensors p eff (at to ^ 0) and C eff are not sign definite. 
This is not at all surprising since ([3]) no longer describes a positive energy for the homog- 
enized dispersive medium. Dispersive effective tensors defined by (|lip correspond to the 
non-local properties of the homogenized medium in the space-time domain where C eff , p eff 
and S cff depend on £ = x — x' and r = t — if . Note that the effective tensors are certainly 
not periodic in k (even if taken on the periodic Bloch branches wb (k) = wb (k + g)), and 
hence they are not periodic in £, as well as non-stationary. It is therefore clear that the 
initial periodic medium with non-zero coupling S (x) in ([2]) cannot be seen as resulting 
from PWE homogenization of another periodic material, say with a much finer scale. At 
the same time, by proceeding from ([2]) and arriving at (|10p with the same entries of the 
coupling terms S and — S + , the PWE homogenization confirms the closure of the Willis 
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elastodynamics model. In this regard, note that our notation S + defined as = 

is equivalent to of [15j (see also the footnote below). 

Consider in more detail the conventional situation where the initial periodic medium 
represents a classically elastic solid with zero coupling S (x) = and real-valued density 
p (x) and stiffness C (x) . It is evident from (|lip with S = that p eff and C eff are even 
functions and S eff an odd function of oj and hence of r. Dependence on k satisfies wb (k) = 
ujb (— k) and 

p cS (oj, k) = p cff * (oj, -k) , C cff (oj, k) = C eff * (oj, -k) , S cff (oj, k) = -S eff * (oj, -k) . (44) 

By (|44|) . C cff , p cff are real and S is pure imaginary at k = and any oj. It also follows 
from (|44p that the C eff (£,r) and p cff (£,t) are real and S eff (£,t) is pure imaginary in the 
space-time domairfl. 

Let p (x) and C (x) be even functions of x. This condition on its own implies that p cS 
and C eff are even functions and S eff is an odd function of k and hence of £. Combination 
with (|44p implies that a classically elastic solid with a periodically even profile is charac- 
terized by p cS and C eff that are real and even (in each variable) functions of either set of 
variables (u, k) and (£,t) , while S eff is a real valued odd function of (oj, k) which vanishes 
at k = and it is a pure imaginary odd function of (£,t). 

The case of an even profile is one example of possible symmetry. The general case may 
be outlined as follows. Let the orthogonal transformation R (R" 1 = R T ) be an element of 
the symmetry group of the initial material, such that (x) = Ri m Rj n RkpRiqC m npq ( x ) 
Vx and also p (x) = p (Rx) and C (x) = C (Rx). Then oj (k) = oj (Rk) and 

Pif (w, k) = R im R jn p c £ n (oj, Rk) , 
Cfjkl ( w ' k ) = R im Rj n Rk p Ri q C^ npq (oj, Rk) , (45) 
Sfjk ( w ' k) = RimRjnRkpS^p (oj, Rk) . 

Consequently, the effective tensors p eS , C eff and S eff are invariant to R at k = 0, and are 
invariant at k ^ to those orthogonal transformations which also satisfy Rk = k. 



5.2 Effective impedance 

Assume the homogenized medium with the effective elastodynamics equations in the Willis 
form (jlOp . The problem of finding the effective response to effective sources © averaged 
over the unit cell of the original periodic medium is now seen as a forced-motion problem: 

{f cff ,7 cfr }(x,t) = {(f),<7)}e i C koe - w *> u cfr (x,t) =(u)e^ k ' x -^. (46) 

This observation makes consistent the use of -S+ in © in place of S 1 " = S T in eqs. (7.3), (7.5)]. 
Equation (J44J 3 implies that the term — S cff in the homogenized equations (|10[) can be interpreted as 
S eff+ (u;, — k) which is the formal adjoint of S considered as a differential operator, i.e. S^, see [131 116] . 
The variety of notation in the literature for the coupling S eff does not make things simpler. In this regard, 
notation in [15] agrees with (|10[) and (|44|l ; notation in [T] eq. (3.26)] is similar to (|10[) although the conjugate 
of Mijk (equivalent to Sijt here) does not appear. 
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In particular, uniform amplitudes (f ) and (7) can now be interpreted as Fourier harmonics 
of arbitrary sources in the homogenized medium. Equations (jlOp with reference to (06]) 
can be recast as 

Z eS (u, k) (u) = (f) + i(wS eff+ - k + C cS ) (7), (47) 
where Z cff is the effective impedance matrix with the components in C 3 

Z$ = Cfakjh - w*,(S$ + St%*) - a, 2 pf (=► Z cff = Z cff+ ) . (48) 

Equation (|4"7|) considered in the absence of sources implies 

Z cff (u B ) = 0, det Z cS (uj, k) = 0. (49) 

Equation (|49fl '? is the dispersion equation for free waves in the homogenized medium, and 
as such may be called the Christoffel- Willis equation due to its analogy to the Christoffel 
equation in a uniform medium. Note that, strictly speaking, Z cff has been introduced 
above as the impedance of the homogenized medium. We now show that it can be equally 
interpreted as the effective impedance for the initial periodic material. The consistency of 
such dual interpretation of Z eff is an essential verification of the model. 

For clarity, consider first the case of a scalar wave problem, where the impedance (|48|) 
is a scalar 

Z cH = k • /i cff k - uk ■ (S cff + S cfr *) - u 2 p cS (50) 

and (|49p reduces to Z eff (cj,k) = 0. Substituting the effective material parameters (|38|) in 
([50]) yields the identity 

Z cS (l>, k) = (Z) - w+GyoW. (51) 

Now recall that (Z) — w + G\ w = by ([33a]), and that, barring certain exceptions 

(see §3.4p . (G)~ 1 (cj,k) = is equivalent to the dispersion equation (J2T]) defining Bloch 
branches in the actual periodic material. Thus we arrive at eq. ([3gD, where Z eS = (G)" 1 
is now an equality, and we observe that Z (w,k) = is the same dispersion equation 
providing both the spectrum of eigenmodes for the homogenized medium and the Bloch 
spectrum (in the main) for the actual periodic material. It is instructive to note from 
(|5ip that the effective impedance Z eS generally differs from (Z) = e + Ze defined by the 
statically averaged material constants (see (|3ip ). This observation highlights the fact that 
the PWE dynamic homogenization of a periodic material must necessarily proceed from a 
forced-wave problem. 

The same conclusion follows for the elastodynamic case. Substituting from (146 p into 
(|48p yields the same expression as obtained for (G)" 1 in (|43|) and therefore verifies that 
Z eff = (G) _1 , now an equality. By (|42p . det(G) -1 = is, apart from exceptional cases, 
equivalent to the Bloch dispersion equation. Comparison with (|49p corroborates the desired 
consistency. Note Z cff ^ (Z) in (|43p . which implies the similar observation as noted above 
for the scalar case. 

The symmetry of the effective material tensors (see £15.ip may lead in a standard fashion 
to quasidiagonal or diagonal structure of Z cff . To be specific, consider an example where 
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k 7^ is parallel to one of the orthogonal translations a,- of a 3D simple cubic lattice of 
spherical inclusions of locally isotropic material in a locally isotropic matrix material. Then, 
using elementary considerations of symmetry, the tensors p eff , C eff and S eff must satisfy 
the tetragonal symmetry (4mm) and hence Z cff is diagonal with two equal components. 
This indicates the polarization of the averaged Bloch eigenmodes defined by (09]) i and also 
shows that the Bloch spectrum uj-q (k) defined by eq. (|49|) -? necessarily contains a set of 
branches that are doubly degenerate for any non-zero k || a^. 

A remark is in order concerning the effective properties for k = 0. Normally eqs. (492) 
and (50) imply that det p eS = (p cS = in the scalar-wave case) at the points k = 0, 
u; (0) ^ of the Bloch branches. It is however noted that Z cff (o;,0) = —uj 2 p cS may not 
hold if uj (0) belongs to exceptional points of divergence of the effective tensors (see §§3.4, 
4). For example, assume a model case of periodic C (x) and constant p = pQ. It is clear 
that Z cff (to, 0) + —u] 2 pqI, since the equality would mean that the Bloch spectrum has no 
open stopbands at k = which is certainly not the case for an arbitrary periodic C (x). 
In fact, spatially averaging the pointwise equilibrium equation 

pou 2 u + diver = (52) 

for k = implies that any Bloch mode with uj ^ must have zero mean, (ub) = 0. In 
conclusion, if p is constant while C (x) is periodic, k = definitely belongs to the excep- 
tional subspace where C eff (cj,0) and/or S cff (u;,0) must diverge. A numerical illustration 
of this is given in £[6j 

5.3 Low frequency limit of the density 

The form of the homogenized effective moduli of (jlip while relatively compact belies their 
intricate dependence on the fixed material properties of the periodic system as well as on 
frequency and wave-vector. Here we touch on one small corner of the parameter space, 
the quasistatic limit, with particular attention to the effective inertia. Assume a classical 
initial medium (with S = and positive p(x)). Taking to — > 0, it is clear from (lllh that to 
leading order in uj the effective moduli reduce to their static form, the effective density is 
the average, and S eff =0(cj). The first correction to the effective density is positive-definite 
tensorial, 

p cff = (p) I + co 2 Y, , pnn P(g)p(-g') G (0) [g, g'] + o(u; 2 ) 

where E g , er \o 6 l 0) [g'g'] d *M ~ s")9Wi = S ^ SS "- 

This illustrates explicitly the departure of the homogenized density from that of classical 
elastodynamics, as expected on the basis of prior results for Willis equations parameters 
for periodic systems [21 [15] . 

6 Numerical illustrations and discussion of the homogeniza- 
tion scheme 

Numerical examples are presented for the ID scalar wave problem. All results reported are 
for the periodic medium based on the the three-layer single cell defined by Table 1. The 



3D dynamic homogenization 



18 




Figure 1: The first six Floquet branches for the ID periodic medium with parameters in 
Table 1. 

PWE effective parameters p cS , p eS and S eS are taken on the Bloch branches wg(fc) which 
reduces them to functions of a single variable, k. The first six branches uiB(k) are shown 
in Fig. [H and the corresponding values of effective parameters are displayed in Fig. [2j 

Table 1: Properties of the three- layer unit cell. Values were chosen to ensure an 
unsymmetric cell with non-commensurate layer thicknesses h. 



layer 


A* 


P 


h 


1 


1 


1 


0.37 


2 


7 


2 


0.313 


3 


1/3 


0.5 


0.317 



Setting apart the obvious case of the origin of the first branch which starts with S = 
and p cS = (p) at k = (see §5.3|) . several other features are noteworthy in Fig. [2j In 
accordance with (|44p . ReS" 35 vanishes at k = on all branches, while ImS" 25 starts at k = 
with non-zero values that alternate in sign beginning with the positive value on the second 
branch. The effective density /? cff vanishes at k = on all branches above the first one, 
as expected from the form of the dispersion relation Z (w, k) = 0, which becomes, using 

k 2 p cS - uk(S cS + S cS *) - uj 2 P cS = 0. (54) 

As k grows, p cS (k) exhibits both positive and negative values, following the alternating 
trend shown by ImS"^. The magnitude of p eS remains within the range of the constituent 
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Figure 2: Effective material parameters calculated for the ID periodic medium (Table 1) 
by the PWE method, eq. (|38[) . The scale is altered on some curves for the sake of legibility. 
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material densities, with the exception of the fifth branch on which it is negative with 
max for a range of k where /? m ax is the largest density in the periodic medium 
(Pmax = 2). The calculated PWE /i eff is always non-negative, however, its magnitude is 
exceedingly large near k = for higher order modes, e.g. the fifth and sixth. This can be 
partly understood from the dispersion relation (|54|) which permits n to become large in 



magnitude as k — > 0. The enormous values on some branches, e.g. fj, > 10, 000 in Fig. 2(e) 



may indicate proximity to an exceptional point related to a small value of det Z\ (see £13.41) . 

The magnitude of S eS is intermediate between that of the density and the stiffness. It is 
worth noting that ([3]) evaluated for the ID scalar Willis model (E — > |(/i eff /c 2 + p cS ui 2 ) \u\ 2 ) 
displays negative values on the third and fifth branches. The proper form of the energy 
density for the dispersive effective medium is a delicate topic [18], best left for a separate 
discussion. 

Figure [3] illustrates the emergence of singular effective parameters at k = 0, w(0) ^ 
when the density is uniform, as predicted in §5.2i Although the curves shown are only for 
the second branch, the same dependence is observed in all higher branches. 

An alternative route to homogenization of lD-periodically stratified elastic body was 
proposed in [2]. The basic idea is to compare the logarithm of the monodromy matrix(MM) 
M (y + T, y) , which is the propagator over one period T along the stratification direction 
Y, with the matrix of coefficients of the system of elastodynamics equations for a uniform 
medium. In this way one can identify the dynamic effective material parameters at any 
frequency wave-vector on a Bloch branch as long as the logarithm is well defined. Unlike 
the PWE, the MM method does not require introducing force terms and is not based on 
averaging. It is worth emphasizing that the PWE and MM approaches are not two 'techni- 
cally' dissimilar derivations of the same result but imply different models of a homogenized 
medium. In this light, the more significant is the result that the MM method also nec- 
essarily leads to the effective elasticity of Willis form (classical elasticity will not suffice) 
[2]. Note that the MM homogenization involves the choice of the reference point y = yo 
at which the initial data for the monodromy matrix M (yo + T, yo) is prescribed. The MM 
homogenized media associated with different yo have the same eigenfrequency spectrum 
oj (k) which by definition coincides with the Bloch branches of the actual periodic medium. 
At the same time, the values of MM effective constants depend on the choice of yo- Once 
any yo is fixed, the exact displacement and traction are recovered at periodically distant 
points yo + nT. As a result, the MM homogenization can provide an exact solution of 
boundary-value or reflection problems for a finite structure containing several periods or 
for a periodic half-space in contact with other regions, see [2]. 

As an example, Fig. H] shows the effective parameters of the same material (Table 1) 
which are calculated for one of the Bloch branches (the fifth) by the MM method. The 
calculation uses the formulas (4.2) and (4.7) i_3 of [2] with a standard definition of the 
2x2 propa gato r M ( yo + T, yo) of scalar waves through a trilayer. The two sets of data on 



Figs. |4(a) and 4(b) correspond to two different choices of the reference point yo within a 
period. Both sets reveal the same general features of MM homogenization \19\ [2] such as 
vanishing of p cS and fi cS at k = tt and a pure imaginary value of S eS . The latter stands 
in stark contrast to the generally complex-valued S eS in the PWE homogenization model. 
At the same time, despite the quite disparate natures of the PWE and MM methods, it 



is noteworthy that they do reproduce the enormous values of // eff occurring at small k on 
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Figure 3: The second Floquet branch for the ID periodic medium with parameters as in 
Table 1, except p% = p\ = 1, and p% = (1 + 8)p\ where 5 = 1.01, 1.05, 1.1. The plot shows 
the developing singular nature of p, cS and S cS at k = for constant density. 

the fifth Bloch branch. This may be seen as evidence that such behavior of p, eS is not a 
numerical artifact. 



7 Conclusion 

The principal result of the paper is a set of explicit formulas for calculating the homogenized 
material parameters using PWE expansions of the underlying elastic moduli and density, 
eq. (jlip . The homogenized governing equations retain important physical properties of 
the original periodic system. Thus, they provide (i) the averaged response to averaged 
source(s), and (ii) the eigen-spectrum w(k) for Bloch waves. The effective parameters of 
eq. pip are convergent for virtually all frequency and wave-vector combinations, including 
{w, k} pairs on the Bloch wave branches. The exceptional cases yield singular values of the 
effective parameters and correspond to zero-mean Bloch waves for the scalar wave system. 

Some comments are in order on the PWE homogenization method and its outcome. A 
nonzero strain source is essential for the derivation of the effective parameters, although the 
strain source may be put to zero after the homogenization is complete. Inclusion of both 
body force and strain source is consistent with the findings of [16] that homogenization 
procedures based on Green's functions can lead to under-determined systems unless a 
full set of sources are included from the outset. The PWE homogenization procedure 
derives equations governing the average of the periodic part of the Bloch wave solution. 
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0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 



k 1% kl% 

(a) yo — 0, at edge of first layer (b) yo — hi, at second layer 

Figure 4: Effective material parameters of the ID periodic medium (Table 1) calculated 
by the MM method for the fifth Bloch branch. The monodromy matrix M(yo + T, yo) m 
(a) and (b) is referred to two different locations of the initial point yo within the unit cell. 



The PWE scheme is not the only approach for defining an effective medium for periodic 
system. Another approach which applies for media with ID periodicity is the monodromy 
matrix (MM) homogenization of [2]. The effective governing equations are again of Willis 
form with explicit expressions for the effective parameters that are distinct from the PWE 
parameters, as the numerical example in ^demonstrates. The fact that different sets of 
effective Willis equations can be obtained for the same medium points to a non-uniqueness 
in the definition of dynamic homogenization. The MM scheme was originally proposed as 
a method for layered structures that faithfully retains the exact phase relationship over 
multiple periods, although at only one point within the period [2]. The PWE scheme 
provides predictions for averaged quantities and does not yield pointwise exact solutions. 

We have shown that the homogenization of a classically elastic medium (S = 0) yields 
a Willis effective medium with either complex- valued S eff (u;,k) (satisfying (|44p ) or with 

pure imaginary S cff (o;) = — S cff (— u) ($6|). The results indicate the closure property for the 
Willis model with strictly complex S under homogenization. 

The present formulation, in starting with a heterogeneous set of Willis elastodynamic 
equations, and resulting after homogenization in a set of homogeneous constitutive equa- 
tions of the same form, confirms the completeness of the Willis model. This type of closure 
complements the fact that the Willis constitutive equations are closed under spatial trans- 
formation |20] , Finally, we note that the present results only partially open the door on 
our understanding of the effective dynamic properties of phononic crystals. Much remains 
to be investigated regarding analytical aspects, numerical implementation and applications 
of these properties. 
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